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Abstract 

We present an approach to simulating quantum computation based on a classical model that 
directly imitates discrete quantum systems. Qubits are represented as harmonic functions in a 2D 
vector space. Multiplication of qubit representations of different frequencies results in exponential 
growth of the state space similar to the tensor-product composition of qubit spaces in quantum 
mechanics. Individual qubits remain accessible in a composite system, which is represented as a 
complex function of a single variable, though entanglement imposes a demand on resources that 
scales exponentially with the number of entangled qubits. We carry out a simulation of Shor's 
algorithm and discuss a simpler implementation in this classical model. 

PACS numbers: 03.65. Sq, 03.65. Ta, 03.67.Lx 
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I. INTRODUCTION 



Quantum computation promises exponential speed-up over classical computation for cer- 
tain problems, such as period finding and quantum simulation. Traditional classical simu- 
lation of a composite quantum system requires updating each of the 2 N amplitudes charac- 
terizing the state of N qubits, according to a Hamiltonian made up of 2 N x 2 N elements. 
The exponential growth of the state space with iV imposes a severe burden on resources for 
this type of simulation. Finding classes of quantum computations that can be simulated 
efficiently is an active field of research 

In quantum mechanics, the qubits that comprise a composite system remain accessible, 
and it is through interactions with individual and pairs of qubits that computation is im- 
plemented. For example, a single-qubit transformation affects all 2 N computational basis 
states of an iV-qubit system. It therefore seems that one of the features of quantum systems 
that enables more efficient computing is the ability to harness the degrees of freedom of a 
2^0 vector space by interacting with only iV qubits. 

Here we present a classical model of discrete quantum systems that is based on repre- 
senting individual qubits and transformations applied to them. This enables us to address 
specific qubits in a composite system, the state of which is represented as a complex function 
of a single variable, and thereby directly replicate the steps in an algorithm as they would be 
implemented in a quantum system. While the resources required to carry out a computation 
exactly are comparable to other methods, this model may be compatible with new approx- 
imations that would enable simulating more qubits than currently is feasible. Furthermore, 
the approach of building a classical model based on imitating quantum systems could offer 
an opportunity to gain insight into the difference in computational power of classical and 
quantum architectures. 

II. REVIEW OF QUANTUM COMPUTATION 

A qubit can be realized with any two-state quantum system that can be prepared in a 
general superposition of basis states of a 2D, complex vector space. For multiple uncoupled 
qubits, the state of the composite system is given by the tensor product of the individual 
states, and the number of computational basis states grows exponentially with the num- 
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ber of qubits. For example, the state of an iV-qubit system, with each qubit in an equal 
superposition of computational basis states |0) and |1), is given by the state vector, 

N 



m = (^=) (io) + ii)) 1 ®(io) + |i)) 
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1 ^ 
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(|0> + |1» 8 ®...® (|0) + |1)) JV 

(|00...00) + |00...01) + 

|00...10) + ... + 111...11)). (1) 



Application of a unitary operation to a qubit that is part of a composite system, which 
constitutes a step in an algorithm, affects all states in the superposition simultaneously, 
illustrating the massive parallelism inherent in quantum computation. 

Any quantum algorithm can be approximated arbitrarily closely using just single qubit 
operations and a generic two-qubit interaction, such as the controlled-NOT (CNOT) gate 3|. 
The effects of these operations can be visualized as rotations and inversions of the 2 N - 
dimensional quantum-computer state vector. The state in Eq. (JTJ), constructed as a product 
of individual qubit states, is a special case. Almost all of the states the system can occupy 
in its vector space will be non-separable, implying that entanglement is required for general 
computation. 

Qubits and operations on them are subject to perturbations from the environment and 
experimental imperfections. In general, it is believed that errors due to decoherence grow 



exponentially with the number of qubits in a system [3|. Realizable quantum computation 
relies on the ability to diminish the effects of these errors-quantum error correction and fault- 
tolerant quantum computation exploit entanglement and the discrete nature of quantum 
systems to make this possible. 

This brief introduction to the fundamental elements of quantum computation emphasizes 
the role played by the mathematical structure of the single- and multiple-qubit vector spaces. 
Quantum systems require this mathematical description, and our classical model is developed 
according to this description by building into it the same state-space structure. The result is 
a new method of simulation and an architecture that may offer insight into the fundamental 
advantages of quantum systems for computing. 
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III. INDIVIDUAL QUBITS 



For the representation of a qubit we use a harmonic function with frequency u. Orthog- 
onal functions sin(u;i) and cos(u;t) serve as convenient basis states and span a 2D vector 
space. We assign these basis functions the role of the computational basis states of a qubit, 

|0) sin(W) 

|1) coa(ut). (2) 

Applying a general unitary transformation U to a qubit i/j, 

ip(t) = asin(a;t) + (3cos(uit), (3) 

requires isolating the coefficients a and j3; each coefficient can then be multiplied by the 
corresponding transformed basis function, yielding the transformed state: 

U[ifj(t)] =aU[wn(vt)] + f3U[cos(u;t)]. (4) 

Orthogonality of the basis functions makes it straightforward to isolate the coefficients by 
taking the inner product of the corresponding function with ip{t): 



LU 

a = — 

7T 



sm(ujt')ip(t') dt', 



o 



,2jt/w 

= - / cos(ujt ; )^(t') dt'. (5) 

JO 

A general transformation is illustrated in Fig. Ufa). 

The modulus squared of the coefficient, or amplitude, of a basis function gives the cor- 
responding "measurement probability". The function representing a qubit can be replaced 
with one or the other basis function according to these probabilities in order to represent 
the measurement process. The ability to determine both measurement probabilities for a 
given state eliminates the need to introduce and carry along normalization coefficients-the 
relative probabilities can be determined at the time of measurement. 

IV. COMPOSITE SYSTEMS 
A. 2 N ~D Vector Space 

This model can be extended to composite systems by using a different frequency for the 
basis functions for each two-state system, creating a new 2D vector space for each qubit. 




FIG. 1: Schematic showing a linear transformation on single and multi-qubit states, (a) For a 
single qubit, orthogonality of the basis functions enables isolation of the coefficients, which can 
then be multiplied by the transformed basis functions, (b) For a composite system, the generator 
G^ n \x,t) enables the functional form of F^ n \t) to be transferred to a different variable. This 
procedure is analogous to addressing a qubit in a quantum system. 

The mapping of quantum states to functions for composite systems becomes 

|0) n sin(u; n t) 

|l) n <^=^ cos(u n t), n = 0, 1, N - 1, (6) 

where n refers to the nth qubit and iV is the total number in the system. 

If N single-qubit functions in equal superpositions are multiplied, the result is a linear 
combination of 2 N different products, 

^N{t) = ( sm(iJit) + cos(uit)) sm(uj 2 t) + cos(uj2t)) ... [sin(ij j^t) + cos(cu j^t)) (7) 
= ( sin(a>it) sin(a;2^) • • • sin(u;Art)) + ( sin(u;i£) sin^i) • • • cos(u;jv£)) + . . . 
+ ( cos(a;it) cos^t) . . . cos(u>]yt)) . 

This is analogous to the tensor-product state for a composite quantum system in Eq. ([I]), 
with the mapping 

\bib 2 b 3 ...b N ) h 1 (u 1 t)h 2 (uj2t)h 3 (u 3 t) . . . h N (u N t) 

= H N;j (t). (8) 
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Here, b n is the binary value representing the state of the nth qubit of a quantum system, h n 
is the basis function (sine or cosine) representing the nth qubit in our model, and H^j{t) is 
the jth of the 2 N combinations of products of h n . 

The functions H^(t) that naturally arise when representing composite systems look like 
iV-qubit computational basis states, and we would like to determine whether they too span a 
state space that grows exponentially with qubit number. This will be the case if all of the 2 N 
functions are orthogonal. It is easy to see by expanding the products of harmonic functions 
in terms of sum and difference frequencies that, for N qubits, the H^{t) are comprised of 
2 N ~ l Fourier frequencies Vti = J2n=i a i,n^n, where is 1 or —1. The Vti are all multiples 
of a fundamenta. frequency given by the greatest common divisor of the qubit frequencies, 
^fund = gcd(o;i, a>2, • • • , oon) pi- The orthogonality of the individual Fourier components 
can be shown to lead to orthogonality of the H^it) for few qubits; the three-qubit case is 
demonstrated in Appendix [A] 



1. Orthogonality of the Hj^(t) 



To demonstrate orthogonality in the general case when the 2 N 1 Fourier frequencies 

are unique, we consider the inner product between two iV-qubit functions H^j(t) = 
h jtl (uit) . . . h jtN {uj N t) and H N>k (t) = h^ujxt) . . . h k)N (u N t), 

e7r/f2 fund 



H NJ (t) ■ H N ,k(t) oc 



%i {ujit) hk, i (uxt) h jt2 (uj 2 t) hk,2 (urf) 

h jtN (u N t)h ktN (u N t)\dt, 



(9) 

where factors for a given qubit have been grouped together Each product 

hj >n {u) n t)h^ n {u) n t) can be written as a function of frequency 2u n , as either | sin(2u; n £) or 
|(1 ± cos(2u; n t)). The integrand in Eq. (jUJ) can then be seen to consist of three types of 
terms. First, there will always be a term that is a product of a function for each qubit, 
hi(2uit) . . . hN{2uj^t). This can be written as a sum of Fourier components with frequen- 
cies that are twice those of the H N {t) and are, therefore, also multiples of VL iunA . Because 
integration of a harmonic function over an integral multiple of its period yields zero, this 
term's contribution to the inner product H^j{t) ■ H^ k {t) vanishes. 

Second, there can be terms that only include factors for some qubits, such as 
h\{2ujit) . . . h n _i{2u) n -ib)h n+ i{2uj n+ it) . . . h N (2u>Nt), which arises when h 
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Fourier frequencies for these terms are also multiples of ilf un( }-the fundamental frequency 
is the greatest common divisor of the set of all qubit frequencies, and necessarily divides 
any subset of them. These terms therefore also vanish when integrated over an interval of 
27r/f2 fund . 

Finally, the integrand in Eq. can have a term of unity (times 1/2 N ). This only occurs 
when all qubit functions are the same for H^ } j(t) and ifjv^t), i.e. HNjif) = H]y^{t). 

2. Redundant Frequencies 

This argument is only valid if all of the Fourier frequencies fli are unique. If there are 
at least two combinations of qubit frequencies that give the same Fourier frequency, we can 
write ^2n=i a n,i UJ n = J2n=i a n,m l ^n, for &i = ^m- Terms that enter this equation with the 
same sign for each Fourier frequency cancel, and the remaining terms give u a + Ub + . . . = 
uo a + ujp + . . ., where the frequencies have been arranged so that all signs are positive. When 
considering all inner products Hpjjit) -HN,k{t), all possible combinations of harmonic factors 
in the integrand in Eq. (Q will arise; for some inner product, there will be a term in the 
integrand like h{2uj a t)h{2ujbt) . . . h(2u a t)h(2oupt) . . ., with Fourier frequencies that include 
Q = 2{uj a + u h + . . . — u a — up — . . .) =0. For the case(s) in which this frequency is the 
argument of a cosine, the constant term results in a nonzero integral, and the different H^{t) 
in this case are not orthogonal. 

Therefore, for sets of qubit frequencies {u> n } that result in 2 iV_1 unique Fourier frequencies, 
the functions H^(t) that naturally arise when representing composite systems are orthogonal 
and span a 2 7V -dimensional space. Unique Fourier frequencies can be ensured by using a 
qubit-frequency definition such as u n = u/2 n ~ 1 . 

B. Linear Operations 

In quantum computation, single-qubit operations along with a generic two-qubit interac- 
tion, such as a CNOT gate, are universal. The general approach to implementing operations 
in our model of composite quantum systems is the same as for a single qubit-we need to iso- 
late the factor multiplying each basis function (sm(uj n t) or cos(u; n i)) for a particular qubit. 
In this case, these factors will be expressions involving other qubit basis functions. Once 
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isolated, they can be multiplied by the transformed basis functions and these recombined to 
generate the transformed composite function. 

Consider a general ^^(t), similar to Eq. (J7J) but with arbitrary coefficients for the if/v(t). 
We can write \I/jv as a sum of two parts, one with terms that include cos(u n t) and one with 
terms that include sin(a; n t): 



where h is cosine or sine, and a^, bk are coefficients for the cos(uj n t), sin(u; n t) terms. FP \t) 
and are the functions that we need to be able to isolate to apply a linear transfor- 

mation to qubit n; determining these functions can be considered to be "addressing qubit 
n." 

For small N, a procedure similar to the one for a single qubit can be adapted. Multi- 
plication of by the relevant basis function for the qubit leads to a different spectrum 
for terms containing that basis function. These different frequencies could be selected from 
the terms containing the orthogonal basis function, enabling the qubit to be addressed. As 
the qubit number grows, however, the number and density of frequencies grow dramatically, 
making this process unfeasible. 

A more general procedure can be used which determines (-Fj ri ' ) ) exactly using the 
orthogonality of the Hiy(t). An inner product can be imposed between and a projector 
that forces all of the terms with one basis function for qubit n to vanish while preserving 
the others. The construction of the projector for a given system is straightforward. 

A linear combination of all of the H^{t) for a system of N qubits can be generated by 
putting each qubit into an equal superposition of basis functions as in Eq. (JTj). We define 
a similar function, the generator for qubit n, as the product of equal superpositions of 
computational basis states for all qubits in the system except qubit n: 

Gjq (t) = ( sin(co>it) + cosfyit)) ( sin(co> 2 £) + cos^i)) • • • ( sin(u;„,_ii) + cos(a;„_it)) 




Fi n \t) cos(u n t) + Fj; n \t) Bin(a; n *). 



(10) 



The H k N {t) are products of harmonic functions 



hi{uit)h 2 {uj2t) . . . h n -\{u) n -ib)h n+x {u) n+ xb) . . . h N (u N t), 



(sin(u; n+ it) + cos(u; n+ it)) . . . (sin(c<jArt) + cos 



(u N t)). 



(11) 



S 



If this is multiplied by cos(a; n t) (sin(u; n t)), the resulting function's inner product with ^jvtX) 
gives the sum of the amplitudes of the terms in Fc n \t) (Fs n \t)Y but the functional form 
is lost. We can salvage the functional dependence by transferring it to a second variable 
introduced to exactly replicate the dependence on t: 

G%\x,t) = 

(sin(a;it) sin(u;ia;) + cos(co>ii) cos(a;i:r)) (sin^t) sin^x) + cos(co>2^) cos^^)) • • • 
( sin(a; n _ 1 t) sin(u; n _ 1 x) + cos(kJ n _it) cos(c<j n _ 1 x)) ( sm(u n+ it) sm(uj n+1 x) 

+ cos(u; n+ it) cos(u; n+ ix)) . . . (sin(c<jjvt) sin(ujNx) + cos(uNt) cos(ujnx)) 
= cos(c<Jit — uj\x) cos(u 2 t — u 2 x) . . . cos(u; n _it — uj n -ix) . . . 

cos(a; n+ it — u n+ ix) . . . cos{u N t — UJnX). (12) 

When this generator is multiplied by cos((x> n t) (sin(u; n t)), we get the projector for 
F^\t) (Fj n) (t)): 

pW(x,t) = cosKt) G%\x,t) 

(pj n) 0M) = sinKt)G^(x,t)). (13) 

Taking the inner product of ^^(t) and Pc H \x,t) (Ps(x,t)), integrated over £, gives us 
F c (n) (x) (Fi n) (x)): 

rir/U {und 

F^(x) = (2 N Q iund /n) / Pj; n \x,t)* N (t) dt 

Jo 

[F^(x) = (2 N Q iund /7c) J Pl n \x,t)* N {t)dt). (14) 

This is an exact technique for addressing a qubit that is part of a composite system. The 
rest of the procedure for applying a transformation follows as in the single-qubit case and 
is illustrated in Fig. [0(b). 

Iteration of this technique of addressing qubits allows for multiple-qubit gates. For ex- 
ample, for a controlled-NOT gate with qubit n\ as the control and n 2 as the target, the 
gate would begin with determination of Fc ni \ which would then take on the role of the 
function \1/ for determining Fc™ 2 ^ and i* 1 ]™ 2 ^. pj™ 1 ^ would be reconstructed after inverting the 
basis functions for n 2 , giving Fc n2 \t) sin(u; 2 ;£) + Fj n2 \t) cos(oo 2 t). Finally, the transformed 
state would be generated as (^Fc n2 \t) sin(co> 2 t) + Fg n2 \t) cos(c<j 2 t) j cos (u> it) + Fg ni \t) sin(c<j 1 t). 
Gates involving more than two qubits can be implemented by further iteration. 
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The measurement probability for a basis function is determined by J^ ntund _p( n )*(£) . 
F {n) (t) dt [6|. Due to the orthogonality of the H N (t), all cross terms in the inner product 
vanish and the result is the sum of the moduli squared of the amplitudes of all of the terms 
containing the corresponding basis function. 

C. Scaling of Required Resources 

The state of a general composite system can be represented as 

uj 2 , u N J?p(uJ Ne+1 )ip(uJ Ne+ 2) • • • , (15) 

where if) represents an individual, unentangled qubit, and ^ characterizes N e entangled 
qubits. While unentangled qubits can be stored and processed individually and with little 
overhead,the resources required to exactly represent entangled qubits scale exponentially 
with N e [7|. The qubit frequencies and the maximum Fourier frequency can be kept finite, 
but the fundamental frequency decreases at least exponentially with number. The interval 
over which \l/ needs to be defined is given by the integration interval required for addressing 
a qubit. In Eq. ffT4j) . the integration limit of 7r/f2f un d yields exact values for the F^; 
coupled with the necessary resolution imposed by the highest Fourier frequency, on order of 
^max/^fund points are required to define 

We can consider the impact on addressing a qubit-both for a unitary transformation 
and for measurement-of simply truncating all of the functions that arise in a calculation. 
For unitary transformations, we determine the accuracy of F^ n \r,t) and F s {n) (r,t), the 
functions determined by integrating Eq. (TH|) to r rather than to 7r/f2 fund , by comparing 
Fc n \r, t) cos(u n t) + Fs n \r,t) sm(uj n t) to ^(t). We define a parameter 8^(t) to represent 
the error in F^ n \r,t): 

+ P( n) (x,t)^(t) dt) sm(u n x) ) 2 dx, (16) 

where ^ is ^ normalized over the interval to r, and N(r) is the normalization constant 
for F c (n) (r, t) cos(uj n t) + F s (n) (r, t) sm(u n t) over the same interval. 



In Fig. [2]^a) we plot S(r) for the state ty(t) = Yl n =i cos ( c< - , n^) + rin=i ^(^n*)) with 



OJr. 



cu/2 n 1 and N e — 5 through 9. We evaluate 5 for the case of addressing the first qubit. The 
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t (units of 27i/co) t (units of 2ji/co) 



FIG. 2: (Color online.) (a) Semi-log plot of error 5 versus integration time r when addressing 
qubit n = 1 for the state discussed in the text. The equal spacing between the steep parts of the 
curves indicates an exponential scaling of integration interval in order to avoid significant errors, 
(b) Plot of r versus integration time r, for the same state and addressing the same qubit. For both 
plots, curves from left to right correspond to N e = 5, 6, 7, 8 and 9 qubits, and each curve extends 
to 7r/r2f un d for the corresponding N e . 

curves show that determination of and Fs n \ which is exact for an integration limit of 
7r/fif un d, abruptly becomes less precise as the integration interval is reduced. Any calculation 
involving many gates will likely require a value for 5 on the steep part of the curve, imposing 
an integration interval that scales exponentially with N e . 

To assess the effect of truncation on measurement probabilities, we truncate the integrals 
used to determine Fc and F^ as above, and then integrate the square of each. We look 
at the ratio of the truncated probability to measure sine versus cosine as a function of 
integration time: 

rW(r) = | (J P( n \x,t)^(t) dt) 2 dxj J (jf P c (n) (MW) d*) 2 dx. (17) 

This ratio is plotted in Fig. [2(b) for the same state used above, for which the actual ratio 
is one. The graph shows that reaching the actual ratio of P s to P c also scales exponentially 
with iV e . However, precise values for measurement probabilities are often not needed and 
useful qualitative information can be obtained for integration intervals that are shorter than 
7r/f2f un d. For instance, for this example the curves indicate that for integration limits beyond 
7r/ (4fi fund ), the ratio r is about 0.25, of the same order as the actual value. 
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V. SIMULATION EXAMPLE FACTORING 



A. The Quantum Fourier Transform and Shor's Algorithm 

The most celebrated quantum algorithm is Shor's method of finding the prime factors of 
a number N. The problem of factorization can be related to the problem of determining the 
period p of the function f(x) = a x (mod N), for an integer a that is co-prime with N; if p is 
even, either (a P//2 + 1) or [dpi 2 — 1) will have a common factor with N [8|. Shor's algorithm 
relies on application of the quantum Fourier transform (QFT) to f(x) to efficiently determine 
the period. 

The qubits involved in implementing this algorithm are divided into two registers, the 
states of which are treated as integers according to the states of the associated qubits. 
The first step in the procedure is to prepare the first register in a superposition of all 
computational basis states, 

2^i_i 

by applying a Hadamard transform to each of the N\ qubits in the register |9j. A gate U 
applied to both registers produces the entangled state 

2*1-1 

\x) ® K(mod N)). (18) 

2=0 

If the second register is measured, it collapses to |a x °(mod N)) for some xq. The first register 
ends up in a superposition of all states \x') for which a x '(mod N) = a x °(mod N), leaving the 
system in the state 

(\x ) + \x +p) + \x + 2p) + . . .) ® |a x °(mod N)). (19) 

Application of the QFT to the first register imposes interference that results in a superposi- 
tion of states \x) that are close to integral multiples of the inverse period, x ^ x K = K,(2 Nl /p), 
for integer k. Measurement yields an integer close to one of the x K , and after several itera- 
tions the period p can be determined. 

The minimum number of qubits required for each register is N± = log 2 N 2 and N 2 = log 2 N. 
This ensures that the second register is large enough to represent p, which satisfies p < N, 
and that the first register is large enough to give a unique value for p from the QFT (see 
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FIG. 3: (a) Illustration of Shor's factoring algorithm for N = 21. The first gate shown represents 
an application of the Hadamard transform to each of the nine qubits in the first register. The 
solid vertical lines that terminate the circuit flow represent measurement; the dashed vertical lines 
indicate points at which ^ is plotted in Fig. UJ (b) Circuit for implementing the quantum Fourier 
transform shown in (a). The unitary gate Rd is a rotation by the angle 7r/2 d , and H is a Hadamard 
transform. The QFT reverses the order of the qubits. 

reference We apply our model to the factorization of N = 21 using a = 2, the first 
integer co-prime with 21. This requires 14 qubits, nine for the first register and five for the 
second; the alogorithm for our example is illustrated in Fig. [3l The qubits are labeled 1 
through 14, with the convention that qubit 1 (14) is the most (least) significant qubit for 
the first (second) register. Qubit n is represented using the frequency u n = u/2 n ~ 1 , and 
the function \I/ representing the state of the system is defined with a resolution of 1/u over 
an interval of 27r/f2f unc j = 27r(2 nmax ~ 1 ) [10]. In Fig. [4] we show the function representing the 
state of the system at various stages in the calculation, as denoted by the dashed lines in 
Fig. da). 

We implement U by generating the output state, shown in Fig. 11(b), by summing all 
functions representing \x) la 1 (mod N)) over x = to x = 2 14 — 1 11]. From here, we measure 
qubits 10 through 14 by comparing Fs^* ■ Fs to Fc* • Fc and applying the rule that 
the outcome of a measurement is the state with the higher measurement probability, or is 
chosen at random if the probabilities are equal. This leaves the second register in the state 
1 10000) = 1 16), and the first register in a linear combination of all \x') between and 2 9 — 1 
for which 2 X ' '(mod 21) = 16; the function representing the first register at this stage is shown 
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FIG. 4: (Color online.) Function ^ representing the state of the system at various stages in the 
simulation, plotted over one period. The different durations shown are due to the different number 
of qubits at those points in the calculation. Only the real part of the function is plotted in (d). 

in FigH](c). Qubits 10 through 14 are removed from the calculation and the QFT is applied 
to the remaining nine qubits. 

The results of applying the QFT and subsequent measurement of qubits 1-9 are shown 
in Fig. Ufa), which displays the measurement probability for the state \x). The peaks at 
multiples of 85| indicate a period of p = 6 for the function 2 x (mod 21), which in turn yields 
either of the prime factors of 21 from gcd(a p / 2 + 1, N) and gcd(a P//2 — 1, N). Figures Efb)-(d) 
zoom in on the probability distributions near different multiples of 85 1. When k x 85 1 is 
not an integer, the probability is distributed among integers closest to the fractional value, 
as in (b) and (c). 

B. Simplifications in a Classical Model 

There is a dramatic simplification to the factoring algorithm when using a classical model 
for simulation. The need for the QFT stems from the fact that, in a quantum system, 
measurement of the state in Eq. (|T9|) results in one of the states in the superposition, with 
no opportunity to learn the others. Repetition of the algorithm to that point likely results 
in a different xq, preventing the period from being learned in successive iterations. In a 
classical system, the state in Eq. f[T9l can be measured as many times as necessary to 
determine the period p. If we apply this simplified scheme to the state \1/ in Fig. HKc), we 
find equal probability for any x' satisfying 2 z '(mod 21) = 16; the first eight peaks are shown 
in Fig. [5](e). 
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FIG. 5: (Color online.) Results of the QFT applied to the first register, (a) Plot of measurement 
probability for all x between and 511. The six peaks are designated by the integer k = through 
5. (b)-(d) Details of peaks for k = 1, 2 and 3 in (a). The dashed line indicates the exact value of 
k(2 Ni /p) for that peak, (e) Plot of measurement probability versus x' for the case when the QFT 
is not used. Only part of the entire distribution, which extends to x' = 511, is shown. 

In addition to avoiding all of the gates required for the QFT, no imaginary numbers are 
required, and most importantly, we save on qubits. The first register, whose size for Shor's 
algorithm is dictated by the QFT stage, in this case only needs enough qubits to represent 
p, and p < N. This provides a savings on order of log 2 N qubits compared to the quantum 
case. Applying this simplified factoring algorithm to iV = 21 would require only 10 qubits, 
five for each register. 



VI. CONCLUSION 



We have presented a classical, qubit-based model of discrete quantum systems that offers 
a new framework for simulating quantum computations by providing access to individual 
qubits and their interactions. The dimension of the state space for a composite system 
grows exponentially with qubit number, and individual qubits continue to be accessible. 
Application to Shor's algorithm highlights the features of the model in action. The resources 
required for implementation scale exponentially with the number of entangled qubits, yet it 
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is possible to save on qubits in a classical model. 



APPENDIX A: FOURIER DECOMPOSITION OF H N (t) 

The ifjv(i) can be Fourier decomposed by expanding the products of harmonics h n in 
terms of sum and difference frequencies. Here we explicitly show the decomposition for the 
case of 3 qubits. The 8 functions H 3 j = i 8 (t) can be expanded as 

COs(cJit) COs(tU 2 t) COs(cU 3 t) = - ^ COS [(Ui + C0 2 + ^3)*] + cos [(ui + U) 2 — ui 3 )t\ 

+ cos [(ui - cu 2 + co 3 )t] + cos [(cji — u 2 — cu 3 )t] j 
cos(c<j 1 t) cos(w 2 ^) sin(o; 3 t) = - ^ sin [(cui + uo 2 + oo 3 )t\ — sin [(cui + uo 2 — tt^i] 

+ sin [(cji - u 2 + u) 3 )t] - sin — u 2 — u 3 )t] j 
cos(cj 1 t) sm(cu 2 t) cos(dJ 3 t) = - ( sin [(a^ + cj 2 + ^3)^] + sin [(cu 1 + u 2 — ui 3 )tj 



- sin [(ui - 0J 2 + cu 3 )t] - sin [(a;i — u 2 — w 3 )t\ j 
cos(c<j 1 t) sin(cj 2 ^) sin(uj 3 t) = — ^ — cos [(0*1 + u> 2 + ^3)*] + cos [(cui + u 2 — u 3 )t\ 

+ cos [(ui - cu 2 + cu 3 )t] - cos [(ui - co 2 - co 3 )t]^ 
sm(ivit) cos(u 2 t) cos(uj 3 t) = - ^ sin [{oui + u 2 + oo 3 )t\ + sin [{oui + uj 2 — uj 3 )tj 

+ sin [(ui - u 2 + cu 3 )t] + sin [(ui - u 2 - u) 3 )i\ j 
sin(c^it) cos(uj 2 t) sin(cj 3 t) = - ^ — cos + u 2 + cj 3 )t] + cos [(oui + uj 2 — uj 3 )t\ 

- cos [(ui - uj 2 + cu 3 )t] + cos — u; 2 — cu 3 )t]^ 
s\n{uit) sin(c<j 2 t) cos(cj 3 t) = - ^ — cos + uj 2 + ^ 3 )t] — cos + uo 2 — oo 3 )t\ 

+ COS [(uJi - U0 2 + UJ 3 )t\ + COS [((jUi —u 2 — oo 3 )t]^ 

sm(uit) sin(cj 2 t) sin(cj 3 t) = - ^ — sin [(cui + uo 2 + c^ 3 )t] + sin [(cu 1 + u 2 — u 3 )tj 

+ sin [(ui - cu 2 + cu 3 )t] - sin [(ui —u 2 — uj 3 )t] j . (Al) 
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We introduce a more compact notation to represent the 8 Fourier components: 



M 0"! 


1 

= 4 


PO^ \ ({i)t -\- dJn -\- (t)r, )f) 


(a 1 n n n n 


1 

= 4 


POS 1 (fill -1- f'fJn fiJo \t\ 


fn n 1 n n n n 

\ ? 5 ? 7 7 5 5 7 


1 

= 4 


POM \ \(iU /j)n -1- (i1n\t\ 


fn n n i n n n 

V 5 5 7 5 5 5 5 / 


1 

= 4 


v v^l ^2 <3 y ^ y 


n n n 1 n n 

^u, U, U, U, 1, U, U, UJ 


1 

= 4 


o i n [it i _l_ / i _l_ / i It 1 

Sin l^c^i ~r u>2 + ui 3 )t) 


(0,0,0,0,0,1,0,0) 


1 

= 4 


sin ((cji + uo 2 - co 3 )t) 


(0,0,0,0,0,0,1,0) 


1 

= 4 


sin ((c^i - oo 2 + ^3)^) 


(0,0,0,0,0,0,0,1) 


1 

= 4 


sin ((cji - u; 2 - ^3)*) 



For all of these terms, (l a ) ■ (I5) = 5 a b(7r/16f2f un d), where (l a ) corresponds to the Fourier 
component represented by the row vector with a one in the ath place and zeros everywhere 
else. 

The functions H 3 (t) written with this vector notation become 

cos(a;it) cos(u> 2 t) cos(w 3 t) = (1,1,1,1,0,0,0,0) 

cos(u;it) cos(u; 2 t) sin(c^ 3 t) = (0,0,0,0,1,-1,1,-1) 

cos(o;it) sin(cj 2 t) cos(u 3 t) = (0,0,0,0,1,1,-1,-1) 

cos(a;it) sin(a; 2 t) sin(c^3t) = (—1,1,1,-1,0,0,0,0) 

sin(c<j 1 t) cos(cj 2 t) cos(w 3 t) = (0,0,0,0,1,1,1,1) 

sin(a;it) cos(a> 2 t) sin^t) = (—1,1,-1,1,0,0,0,0) 

sin(a;it) sin(cj 2 t) cos(w 3 t) = (—1,-1,1,1,0,0,0,0) 

sin(u;it) sin(u; 2 t) sin^t) = (0,0,0,0,-1,1,1,-1). 

From this it can be seen that all of the H 3 (t) are orthogonal. For example, 

cos(o;it) cos(u> 2 t) cos(cj 3 t) • cos(u;i£) sin(c«j 2 t) sin^t) = 
(1,1, 1,1, 0,0, 0,0) -(-1,1, 1,-1, 0,0, 0,0) = 

-1 + 1 + 1-1 = 0. (A3) 
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In this case and for small numbers of qubits, all of the different inner products can be verified 
to be zero, and Hnj ■ H^j = (7r/4f2f un d). 

We can also see that not all of the functions are orthogonal if there is redundancy in the 
Fourier frequencies. If, for instance, U\ -\- 002 — ^3 = uji — uj 2 + uj 3 , then the second and third 
vectors in Eq. flA2|) are identical (as are the sixth and seventh). In this case, the H 3 (t) can 
only be represented by six orthogonal vectors, not eight. Then, 

cos(a>it) cos^t) cos(u;3t) • cos(o;it) sin^t) sm(uj 3 t) = 
(1,2, l,0,0,0)-(-l, 2,-1,0,0,0) = 

-1 + 4-1 = 2. (A4) 
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